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Multi-time density correlation functions in glass-forming liquids: 
Probing dynamical heterogeneity and its lifetime 

Kang Kim x 'Eland Shinji Saito 1 ^ 

Institute for Molecular Science, Okazaki 444-8585, Japan 
(Dated: 28 June 2010) 

A multi-time extension of a density correlation function is introduced to reveal temporal information about 
dynamical heterogeneity in glass-forming liquids. We utilize a multi-time correlation function that is analogous 
to the higher-order response function analyzed in multidimensional nonlinear spectroscopy. Here, we provide 
comprehensive numerical results of the four-point, three-time density correlation function from longtime 
trajectories generated by molecular dynamics simulations of glass-forming binary soft-sphere mixtures. We 
confirm that the two-dimensional representations in both time and frequency domains are sensitive to the 
dynamical heterogeneity and that it reveals the couplings of correlated motions, which exist over a wide range 
of time scales. The correlated motions detected by the three-time correlation function is divided into mobile 
and immobile contributions that are determined from the particle displacement during the first time interval. 
We show that the peak positions of the correlations are in accord with the information on the non-Gaussian 
parameters of the van-Hove self correlation function. Furthermore, it is demonstrated that the progressive 
changes in the second time interval in the three-time correlation function enable us to analyze how correlations 
in dynamics evolve in time. From this analysis, we evaluated the lifetime of the dynamical heterogeneity and 
its temperature dependence systematically. Our results show that the lifetime of the dynamical heterogeneity 
becomes much slower than the a-relaxation time that is determined from the two-point density correlation 
function when the system is highly supercooled. 



I. INTRODUCTION 

When liquids are supercooled below their melting tem- 
peratures while avoiding crystallizations, they eventually 
undergo a glass transition to become amorphous solids. 
The glass transition is ubiquitous among a wide variety 
of materials. There are many known properties associ- 
ated with the glass transitions^— In particular, when 
the glass transition is approached, various time corre- 
lation functions decay with non-exponential relaxations. 
Moreover, dynamical properties such as the structural 
relaxation time and the viscosity of the system tend to 
diverge, whereas the static structures remain unchanged 
and thus similar to those of normal liquids. Despite a 
large number of theoretical, experimental, and numerical 
studies over the past decades, the understanding of the 
mechanisms behind this drastic slowing down remains 
one of the most challenging problems in condensed mat- 
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To address this problem from the microscopic level, 
various experiments have been employed including nu- 
clear magnetic resonance and optical spectroscopies.— ~— 
These studies have shown that the dynamics do not fol- 
low the "homogeneous" scenario, but instead follow the 
"heterogeneous" scenario in glass-forming liquids^ - — In 
the heterogeneous scenario, the non-exponential relax- 
ation is explained by the superposition of individual par- 
ticle contributions with different relaxation rates. 

Recent molecular dynamics (MD) simulations of model 
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glass-forming liquids^— and experiments performed 
on colloidal dispersions using particle tracking tech- 
niques^— have provided direct evidence that the struc- 
tural relaxation in glassy states occurs heterogeneously, 
i.e., there is a coexistence of mobile and immobile states 
moving within correlated regions. These studies have also 
shown that the sizes of the correlated regions gradually 
grow beyond the microscopic molecular length scale with 
decreasing temperature (increasing the volume fraction 
in the case of colloidal dispersions). Such recent efforts 
have established the concept of "dynamical heterogene- 
ity" (DH), an idea that advocates for a key mechanistic 
role underlying the drastic slowing down of the glass tran- 
sition. Thus, to understand the details of the relaxation 
processes involved in DH, we must systematically char- 
acterize and quantify its spatiotemporal structures. The 
questions we seek to answer, then, include "how large 
are the heterogeneities?" and "how log do they last?" as 
discussed in Ref. [l2|. 

Recently, the determination of the size and length scale 
of the DH has attracted much attention. The correla- 
tions in dynamics can be measured in terms of four-point 
correlation functions and their associated dynamical sus- 
ceptibility. This approach has been successful in extract- 
ing and characterizing the growing length scale with ap- 
proaching the glass transition! 18 i 24 i 25 ' 28 i 35 ~ — The four- 
point dynamical susceptibility has also been investigated 
by the mode-coupling theory 4 ^— and experiments.— ~— 

However, knowledge and measurements relating to the 
time scale and lifetime of the DH are still limited. More- 
over, the temperature dependence of the lifetime remains 
controversial. It has been observed in some numerical 
simulations that the lifetime and characteristic time scale 
of the DH are comparable to the a-rclaxation time, r Q , 
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as determined by the two-point density correlation func- 
tioni 20 ' 23 i 55 On the contrary, other simulations show that 
the lifetime becomes much slower than the r a as tem- 
perature decreases, and indicate that there exist devia- 
tions between the two time scales in glass-forming mod- 

c fc 19.28.56-59 

The aim of the present paper is to investigate the DH 
by numerically calculating the multi-time density corre- 
lation function, which is an elaboration of our previous 
studyi« In this paper, we emphasize the essential con- 
sideration of the multi-time extension of the four-point 
correlation function that can aid in the elucidation of the 
time evolution of the correlated particle motions in the 
DH. We show numerical results of the multi-time corre- 
lation function via the two-dimensional (2D) represen- 
tation analogous to the multidimensional spectroscopy 
techniques. The 2D representations in both the time and 
frequency domains enable us to explore the couplings of 
particle motions in the DH. Furthermore, the multi-time 
correlation function is divided into mobile and immo- 
bile contributions from the single-particle displacement. 
It is demonstrated that this decomposition provides ad- 
ditional information regarding detailed relaxation pro- 
cesses of both mobile and immobile correlated motions in 
the DH. From extensive numerical results of the multi- 
time correlation function, we determine the lifetime of 
the DH and resolve all controversy regarding the tempo- 
ral details of the DH. 

The paper is organized as follows. In Sec. [ill we briefly 
review recent studies that have used the four-point cor- 
relation function and its associated dynamical suscepti- 
bility to characterize the correlation length of the DH 
in glassy systems. Furthermore, we highlight how the 
multi-time extension is a crucial element in the discovery 
of temporal information of the DH. In Sec. IIII1 we briefly 
review our MD simulations and summarize some numer- 
ical results using conventional time correlation functions. 
In Sec. lIVl we present numerical calculations of the multi- 
time correlation function and the time evolution of the 
correlated motions of the DH. We also determine the 
lifetime of the DH and its temperature dependence. In 
Sec. El we summarize our results and give our concluding 
remarks. 



II. MULTI-POINT AND MULTI-TIME CORRELATION 
FUNCTION 

A. Four-point correlation function to measure the 
dynamical correlation length 

As mentioned in the introduction, the concept of the 
DH indicates that the mobility of individual particles 
largely fluctuate in the slow dynamics. Furthermore, par- 
ticles that have similar mobility form cooperative corre- 
lated regions. Conventional analysis that is based on 
the use of two-point density correlation functions, for 
example the intermediate scattering function F(k, t) = 



(p(k, i)p(—k, 0))£i, cannot detect large fluctuations in lo- 
cal mobility because two-point correlation functions aver- 
age over all particles. Here, p(k, t) = YljLi sxp(ik-Tj(t)) 
is the Fourier transform of the density held of the N par- 
ticles in the system. Tj(t) is the jth particle position at 
time t and k is a wave vector with k = \k\. 

To characterize and quantify the correlations of the 
local mobilities, we need to analyze the correlations of 
the fluctuations in the two-point density correlation func- 
tion^^^r^B. This i eads to the following four- 
point correlation function, 

4 4 )(M) = N{SF q (k,t)SF_ q (k,t)}, (1) 

with 

1 N 

SF q (k,t) = -J2e iq - rm (eMik-&r j (t)]-F(k,t)). 

3=1 

(2) 

Here, q is a wave vector with q = \q\. Integrating over the 
volume and setting q — > 0, we obtain the so-called four- 
point dynamical susceptibility, X-iW- If the DH becomes 
dominant in the slow dynamics and if the fluctuations in 
the particle mobilities become large, then the Xi{t) wm 
be able to show the growth of its correlation length, £. 

Recently, the four-point dynamical susceptibility Xi(t) 
has been intensely applied to study physical implementa- 
tions of the DH in various systems that include sheared 
supercooled liquids,— aging in structural glasses^ su- 
percooled water^ 5 - slow dynamics confined in random 
media^ colloidal gelations^ and sheared granular ma- 
terials j££ 

It is remarked that the value of Xift) depends on the 
choice of the ensemble. This ensemble dependence in- 
fluences the estimation of the correlation length £ for 
q -> 

B. Why use a multi-time correlation? 

The four-point correlation function defined by Eq. (JTJ) 
is a one-time correlation function, as is schematically il- 
lustrated in Fig. HJa). In order to quantify the tempo- 
ral details of the DH and its lifetime Thctoro, it is essen- 
tial to analyze how the correlated particle motions decay 
with time. This requires a multi-time extension of the 
four-point correlation function. In practice, by following 
Fig. [Ijb), the four-point correlation for the density field 
p(k,t) can be generalized to the three-time correlation 
function with correlations at times 0, T\, T2, and t 3 given 
by 

F 4 (ki,k 3 ,t 1 ,t 2 ,t 3 ) 

= (p(k3,T 3 )p(-k 3 ,T 2 )p(k 1 ,T 1 )p(-k 1 ,0)). (3) 

This equation takes into account three time intervals, 
ti, £2 j an d t 3 . The definition of the time interval ti is 
ti = Ti — Ti—i, where tq = 0. 
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(a) 




^1 x 2 ^3 

FIG. 1. Schematic illustrations of the time configuration of 
the four-point correlation function; (a) four-point correlation 
function denoting the correlation of fluctuations in the two- 
point correlation function between two times and t, and (b) 
the three-time correlation function with correlations at four 
times 0, ri, T2, and ts. 

We can define the following function as the difference 
between the four-point and three-time correlation func- 
tion, Fi{ki, &3, t\, t 2 , ta), and the product of the two- 
point correlation functions: 

AF(k u k 3 ,h,t 2 ,t 3 ) 

= F 4 (fei,A;3,ti,i2,i3) - F(ka,ta)F(ki,tx). (4) 

This can be regarded as a multi-time extension of Eq. ([1} . 
If the dynamics are homogeneous and if the motions be- 
tween the two intervals t\ and £3 are uncorrelated and 
decoupled, then the three-time correlation function A.F4 
should become zero. On the other hand, if the dynamics 
become heterogeneous, the dichotomy between the mo- 
bile and immobile regions would lead to finite values of 
AF4 because of the correlated motions between the two 
intervals t\ and t 3 . Furthermore, the progressive changes 
in the second time interval t 2 = T3 — t 2 of AF4 enable us 
to investigate how the correlated motions between two 
time intervals t\ and t 3 decay with the waiting time t 2 . 
This can provide the temporal information regarding the 
DH that is relevant in the quest to quantify its lifetime, 

' nctcro^^ 

Some computational studies have already utilized 
multi-time correlations to examine the heterogeneous dy- 
namics ,19120,23,55,56,69-^2 However, in these calculations, 
only limited information has resulted (for example, re- 
sults for ti = ts = r a has been successfully provided). 
This lack of results has caused the aforementioned con- 
troversy regarding the temporal information that is rele- 
vant to the DH. Throughout this paper, we present the 
comprehensive numerical results of a four-point, three- 



time density correlation function without fixing any time 
intervals. This multi-time correlation function is used to 
quantify the lifetime of the DH, ihetero and to determine 
its temperature dependence. 

It is of interest to note that the multi-time correlation 
function can be regarded as an analogue of the nonlin- 
ear response functions of a molecular polarizability and 
dipoles as analyzed using the multidimensional spectro- 
scopies such as 2D Raman spectroscopy and infrared 
(IR) spectroscopy^ - — There exist promising theoretical 
treatments for the multi-time correlation function based 
on the mode-coupling theory 79 ' 80 These techniques have 
now become powerful and standard tools to study con- 
densed phase dynamics. For example, there are often 
used to study the ultrafast dynamics of liquid water 
The utility of these techniques is enabled by the ability 
of the nonlinear response function to reveal details about 
the couplings between motions. This information is not 
available in the one-time linear response function. The 
present study analogously employs the underlying strate- 
gies and concepts of these multidimensional spectroscopy 
techniques to study the heterogeneous dynamics of the 
glass transition. 

It should be remarked that unique experiments have 
been recently proposed to examine heterogeneous dy- 
namics in various chemical systems, which are referred 
to as 2D Fourier imaging correlation spectroscopy 8 - 9 - and 
multiple population period transition spectroscopy 90 ' 91 
These techniques provide information based on four- 
point correlation functions, which are basically the same 
as Eqs.Oandgl 

III. SIMULATION MODEL AND SOME DYNAMICAL 
CONSIDERATIONS 

A. Model 

We carried out MD simulations for a three-dimensional 
binary mixture. Our system consists of Ni = 500 parti- 
cles of component f and N 2 = 500 particles of component 
2. They interact via a soft-core potential 

v ab {r) = e^ff\ (5) 

where 



and a, b G {1,2}. The interaction was truncated at 
r = i<7\. The size and mass ratios were <J\ja 2 = 1/1.2 
and mi/m 2 = 1/2, respectively. The total number den- 
sity was fixed at p = N\ + N 2 /L 3 = O.Sc-f 3 , where the 
system length was L = 10.77cti under periodic boundary 
conditions. In this paper, numerical results will be pre- 
sentcd in terms of reduced units 01, e/ks, t = \Jm\o\jt 
for length, temperature, and time, respectively. The 
velocity Verlet algorithm was used with a time step of 
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FIG. 2. Time correlation functions of glass-forming liquids 
at temperatures T = 0.772, 0.473, 0.352, 0.306, and 0.289 
from left to right; (a) the self-part of the intermediate scat- 
tering function F s (k,t) with k = 2n, (b) the mean square 
displacement (5r 2 (t)}, (c) the non-Gaussian parameter ot2{t), 
(d) the new non-Gaussian parameter *y(t), and (e) the four- 
point dynamical susceptibility x±(k,i) with k = 2tt. Inset of 
(a): a-relaxation time r a as a function of the inverse tem- 
perature 1/T. The dotted line in (b) refers to the diffusion 
asymptote 6Dt at T = 0.289. 



0.005r in the microcanonical ensemble. The states in- 
vestigated here were T = 0.772, 0.473, 0.352, 0.306, and 
0.289. Below, we present and summarize the key numer- 
ical results regarding the dynamic properties of the glassy 
dynamics. Other information regarding this model, in 
particular static properties such as static structure fac- 
tors can be found in some previous worksi 18 ' 21 



B. Intermediate scattering function and mean square 
displacement 

To begin, we examined the density fluctuation in terms 
of the self-part of the intermediate scattering function of 
component 1 particles that is defined as 



F s (k,t) = (^-J^e X p[ik-A rj (0,t)] 



(7) 



where Atvj(0, i) = rj(t) — Vj(0) is the jth particle dis- 
placement vector during the two times and t. The 
behavior of F s (k,t) is often utilized to study the non- 
exponential decay of the structural relaxation as demon- 
strated in Fig. [2ja). Here the wave vector is chosen as 
k = 2tt. This corresponds to the wave vector of the first 
peak of the static structure factor. Also well known is the 
fact that when the temperature decreases, this function 
plateaus during the /3-relaxation regime. In this regime, 
a tagged particle is trapped by its surrounding caged par- 
ticles. Eventually, the tagged particle escapes from the 
cage on a much longer time scale, which is referred to 
as the a-relaxation regime. In this paper, we define the 



a-relaxation time r a as F s (k,r a ) 



with k = 2tt. In 



the inset of Fig.[2][a), the temperature dependence of r a 
is plotted as a function of the inverse of the temperature 
1/T. We observe that the structural a-relaxation time 
is drastically increased and exhibits super- Arrhenius be- 
havior as the temperature is decreased. 

Particle motions are also analyzed through the mean 
square displacement (MSD). We calculated the MSD for 
the particles of component 1, 



1 Nl 

^X>^(o,*)l 5 



(Sr 2 (t)) 



(8) 



and displayed the results in Fig. &b) for various tem- 
perature. As shown in Fig. [2jb), at lower temperatures 
a plateau develops during the intermediate /3-relaxation 
regime, when the cage effect is dominant. Diffusive be- 
havior, (5r 2 (t)) = 6Dt , eventually sets in over the time 
scale of t ~ t„. 



C. Non-Gaussian parameter 

We next employ the non-Gaussian parameter (NGP) 
a 2 (t) defined as 



a 2 (t) 



3(^ 4 (t)) 
5(Sr 2 {t)) 2 



(9) 



j'=i 



with (SrHt)) = {(l/iVOE^IAT^O,;)! 4 ). The a 2 (t) 
reveals how the distribution of the single-particle dis- 
placement, Sr, at time t deviates away from the Gaus- 
sian distribution^ As is well documented^ and shown 
in Fig. [2jc), a^it) begins to grow as the temperature is 
decreased. The growth of a 2 (t) means that the distribu- 
tion of the displacement will have two peaks. These peaks 
indicate the existence of both mobile and immobile par- 
ticles, the main feature of DH. However, it is noted that 
ct2(t) mainly grows in the /3-relaxation regime, when the 
F s (k,t) plateaus (see Fig. [UJa)). In practice, in the time 
scales at the beginning of the /3-relaxation regime, a 2 (i) 
begins to grow, whereas on the time scale of r ai a 2 (t) 
begins to decrease to zero. This is due to the fact that 
a 2 (t) is strongly dominated by the mobile particles which 
move faster than particles with a Gaussian distribution. 
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FIG. 3. Distribution of the logarithm of the single-particle 
displacement P(\og 10 (6r);t) at T = 0.772 (a), 0.352 (b), 
and 0.289 (b). For each temperature, the times shown 
are t — r a (solid curve), 2r a (dashed curve), 4r a (short 
dashed curve), and 10r a (dash-dotted curve). The dotted 
curve in (c) refers to the Gaussian distribution G s (Sr, t) — 
[l/(47rDi) 3/2 ]exp(-<5r 2 /4Di) at t = 10r a , where the diffu- 
sion constant D is evaluated by the asymptote 6Dt of the 
mean squared displacement (see Fig. [2jb)). 



The time tngp during a^CO has a peak thus becomes 
smaller than t q at lower temperatures. It is remarked 
that a similar behavior in ot2(t) has been demonstrated 
using the mode-coupling theory, which incorporates hop- 
ping motions^ 

Instead of using the NGP ct2(t), Flenner and Szamel 
have recently proposed a new non-Gaussian parameter 
(NNGP), 7(t), defined as 



1 



Sr 2 (t) 



-1, 



(10) 



with (1/Sr 2 (t)) = <(l/JVi)Epi lAr^O,*)!- 3 ). The 7 (t) 
strongly weights the immobile particles which have not 
move as far as the Gaussian distribution would predict^ 
Figure Hd) demonstrates 7(f) for various temperatures. 
It is observed that the time tnngp a t which j(t) has a 
peak is of a longer time scale than r Q . This is in contrast 
to the results of the conventional NGP analysis as shown 
in Fig. He). 



D. Four-point dynamical susceptibility 



As mentioned in Sec. Ill A[ the four-point correlation 
function that is defined as the correlation function of 
the fluctuations in the two-point correlation functions 
has become a powerful tool to determine the correlation 
length of the DH. Although there are several definitions 
for x 4 (/c,i), one is given by^S 



Xi(k > t)=N 1 



1 Nl 

w y £SF j (k,0,t) 



where 



<5Fj(fe,0,t) = cos[fc • Arj(0,t)} - F s (k,t), 



[11) 



(12) 



represents the individual fluctuations in the real-part and 
self-part of the intermediate scattering function between 



time and time t. Alternatively, Xi(k,t) can be ex- 
pressed by22- 

X 4(k,t) = iVx[(F s (fe,t) 2 ) - (F s (k,t)) 2 }. (13) 
Here we adopt F s (k,t) as 
1 Nl 

F s (k,t) = —J2cos[k-A rj (0,t)], (14) 
1 j=i 

with F s (k,t) = (F s (k,t)}. The X±{k,t) shows the cor- 
relation of the fluctuation in the two-point correlation 
function F s (k,t). This reveals how the particle motions 
(or trajectories) between times and t are correlated. 
In other words, the amplitude of Xi{k,t) signals the 
total amount of spatial correlations in the particle dis- 
placements within the given time interval t. As seen in 
Fig. He), the Xi(k,t) typically presents non-monotonic 
time behavior. The peak of Xi{k,t) appears on a time 
scale that is comparable to r Q . Note that the ensemble 
dependence of the dynamical susceptibility is not taken 
into account since the microcanonical dynamics is em- 
ployed in our simulations. 



E. Distribution of single-particle displacements 

We end this section with a discussion of the distribu- 
tion of single-particle displacements, as alluded to above. 
Following Flenner and Szamel^ we calculated the dis- 
tribution P(log w (Sr);t) of the logarithm of the particle 
displacements, Sr, at time t. These displacements are 
obtained from the self-part of the van-Hove correlation 
function G s (5r,t) as 



P(log w (Sr);t) = ln(10)4ir6r 3 G s (5r,t). 



(15) 



Figure [3] shows P(\og 10 (5r); t) at various times t for 
T = 0.772, 0.352, and 0.289. The dotted red curve in 
Fig. He) refers to the distribution of the Gaussian pro- 
cess, G s (Sr,t) = [l/{4irDt) 3 / 2 }cx-p(-6r 2 /4:Dt) with the 
diffusion constant D at T = 0.289. It is noted here that 
the Gaussian distribution is independent of time i, and 
the peak height is given by P(\og w (6r);t) ps 2.13.^ At 
the high temperature T = 0.772, there is only one peak 
during time t. This peak has the smallest deviation from 
the Gaussian distribution. On the contrary, at the lowest 
temperature T = 0.289, we clearly see the two distinct 
mobile and immobile peaks. These peaks clearly have 
large deviations from the Gaussian even for longer time 
scales than r a , implying the long-lived DHiSi 



IV. NUMERICAL RESULTS OF THE MULTI-TIME 
CORRELATION FUNCTION 

A. Three-time density correlation function 

Following the time configuration illustrated in Fig.fTJb) 
and Eq. Q, we extend the dynamical susceptibility 
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FIG. 4. 2D representations of the three-time corre- 
lation functions; (a) total AFi(k,t\,t2,tz), (b) the mo- 
bile part A_F™°(fc, ti, t%, tz), and (c) the immobile part 
A_F4 m (fc, ti, t2, ts) at waiting time t% = for various tem- 
peratures T = 0.772, 0.352, and 0.289 from left to right. The 
wave vector k is chosen as k = 2n. 
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FIG. 5. Diagonal parts of the three-time correlation func- 
tions, AF 4 , AF™°, and Af| m at waiting time t 2 = for 
T = 0.772 (a) 0.352 (b), and 0.289 (c). The solid, dashed, 
and short dashed curves correspond to total, mobile, and im- 
mobile parts, respectively. 



that move more (less) than the mean value of the single- 
particle displacement, */ (Sr 2 (ti)}, for the first time in- 
terval t\. During t\, the function SFj(k, 0,ri) selects 
the sub-ensemble of mobile (immobile) contributions in 
the DH, and the total function SFj(k, T2,Ts)SFj(k, 0,Ti) 
contains information to determine how long the mo- 
bile (immobile) particles remain during the waiting time 
t 2 . This information gained from the three-time corre- 
lation function AF4 is related to the joint probability 
P (Sr(t 3 )\5r(ti)) of two successive particle displacements, 
5r(t 3 ) and 8r(t\). This joint probability describes the 
probability of the particle being mobile (immobile) at t\ 
and remaining mobile (immobile) at t 3 after the waiting 
time t 2 . More details will be given elsewhere by analyz- 
ing the multi-time correlation functions of the particle 
displacements^ 



X<i(k,t) defined by Eq. (fTTj) to the three-time density 
correlation function with times 0, n, t 2 , and t 3 . This 
is defined as 

AF 4 (k,h,t 2 ,t 3 ) 

I 1 Nl \ 

= ^E^( fc ' T2 ' T 3)^( fc ' ' T i))/- ( 16 ) 

We note that Eq. (flB)) is regarded as the self-part of 
Eq. (Q| . Moreover, the wave vector is chosen as k = 
k\ = k 3 in our numerical calculations. As discussed 
in Sec. Ill Bl AFn{k,t\,t 2 ,t 3 ) denotes the correlations of 
fluctuations in the two-point correlation function F s (k, t) 
between two time intervals, ti = T\ and £3 = T3 — t 2 . 

Furthermore, the three-time correlation function given 
by Eq. (|16j) can be divided into two parts as 

AF 4 (k,t u t 2 ,t 3 ) 

= AF™°(k, h,t 2 , t 3 ) + AF™(k, h,t 2 , t 3 ), (17) 

where AF™ (AF] m ) represents the mobile (immobile) 
part, arising from the contribution of mobile (immobile) 
particles during the first time interval, t\. Practically, we 
defined the mobile (immobile) particles as those particles 



B. Zero waiting time t 2 — 

We first present the numerical results of the three- 
time density correlation function, AF 4 (k, t±, t 2 , £3), at the 
waiting time t 2 = 0. These are shown in Fig. [4ja) at 
various temperatures, T = 0.772, 0.352, and 0.289. It 
can be seen that the intensity of AF^k, ti,t 2 ,t 3 ) grad- 
ually grows with decreasing the temperature. This in- 
dicates that particles that are mobile (immobile) during 
the first time interval, t\, tend to remain mobile (im- 
mobile) during the subsequent time interval, t 3 . It can 
also be seen that the profile of AF 4 (k, ti,t 2 , t 3 ) is widely 
broadened, suggesting that the motions between various 
time scales, including a- and a-relaxation and a- and (3- 
relaxation, are coupled. Furthermore, the time at which 
AF 4 has its maximum value is approximately given by 
the a-relaxation time, T a . 

To describe the details of AF^k, t\, t 2 , t 3 ) more fully, 
we show the mobile and immobile parts of the system, 
AFl lw (k,h,t 2 ,t 3 ) and AF| m (fc, t u t 2 , i 3 ) in Fig.^b) and 
(c), respectively. These diagonal parts at t\ = t 3 are also 
drawn in Fig. [5] for T = 0.772 (a), 0.352 (b), and 0.289 
(c). It can be seen that the peak of AF4 is composed 
of the two distinct mobile and immobile contributions 
particularly at the lower temperatures. The time scales 
of the two contributions are different, i.e., the peak of 
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FIG. 6. 2D representations of the three-time corre- 
lation functions; (a) total AFi{k,t\,t2,t-i), (b) the mo- 
bile part AF^°{k, ti, is, £3), and (c) the immobile part 
AF| m (fc,ti,i 2 ,t 3 ) at T = 0.289. Waiting times are varied 
for ti = 0, T a , and 3r a from left to right. The wave vector k 
is chosen as k = 2n. 



the mobile part, AF4 110 , appears at t% ~ tngp> while the 
peak of the immobile part, AF£ 10 , is pronounced on the 
time scale of t\ ~ tnngp- These findings are expected as 
per the discussion in Sec. lIII Cl Specifically, the NGP fo- 
cuses on the mobile particles and the NNGP weights the 
immobile contributions of the non-Gaussian distribution 
of the particle displacement (see Fig. GUc)). 



C. Waiting time t 2 dependence and the lifetime of the 
dynamical heterogeneity 



As outlined in Sec. Ill Bl the progressive changes in the 
waiting time t% = T3 — ti of the three-time correlation 
function AF^k, t\, ti, £3) make it possible to investigate 
how the correlated motions decay with time. Figure [5] 
shows the time evolutions of the three-time correlation 
functions, AF4, AF™° , and AF" n at the lowest temper- 
ature of T = 0.289. We also plot the diagonal parts of 
the evolution along t\ = at various t% in Fig. [7J It 
is demonstrated that the correlations gradually decay as 
the waiting time t% increases. The values of the three- 
time correlation functions tend toward zero for ti — s> 00. 
We also find that for larger ti, the peak of AF4 tends 
to shift to t ~ 1000, which is close to tnngp- This peak 
shift is attributed to the fact that immobile particles tend 




C0.02 



(b) 




(c) 




10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 



FIG. 7. Diagonal parts of the three-time correlation func- 
tions; (a) total A.F4, (b) the mobile AF4 10 , and (c) the immo- 
bile A_F| m at T = 0.289. Waiting times are varied at t 2 = 0, 
0.5r a , Ta, 2r a , 4r a , 6r a , and 10r Q from top to bottom. 



(b) 



io° I LJ 1 °° I ^ 
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FIG. 8. 2D representations of relaxation time distributions 
of three-time correlation functions; (a) fhetero, (b) fhetero i and 
(c) fhetero- normalized by r a at T — 0.289. The dotted line 
refers to the iso-line of Thetero = t q in each panel. 



to remain immobile on larger time scales as indicated in 
Fig. [UJc) . Moreover, the presence of correlations on larger 
time scales than the a-relaxation time scale is observed, 
even for ti = 10r Q . It is also clearly seen that the off- 
diagonal parts of AF4, AF" 10 , and AF™ become notice- 
able with increasing ti. These observations imply that 
the relaxation rates of AF 4 , AF™ , and AF™ largely 
depend on which time scale is examined. 

To explore the details of the time scale of the correlated 
motions, we define the relaxation time fh e tero(ii, £3) of 
the AFi(k, ti, ti, t 3 ) as 



AFi{k,t U Thetero, h)/AF A {k, t u 0,t 3 ) = e 



(18) 



for various values of t\ and t 3 . Similarly, the relaxation 
times T^° tcm (t x , t 3 ) and r h ™ cro (ii , t 3 ) are determined from 
AF" 10 and AF] m , respectively. Figure M shows the 2D 
representations of the relaxation times fhetero, f^toro' anc ^ 
f h ™ ero at T = 0.289. We confirm that the relaxation 
time fhetero of the total function AF<± is described by 

.mo 
hctcro 



the summation of the mobile and immobile parts, f. 



*im 
hctcro' 



the distribution of the fhetero has a multiple structure; 
the relaxation time fhetero becomes much larger than the 
time scale, t q , if the time interval t\ or the time interval 
t 3 is examined for larger time scale than t q . On the 
other hand, fhetero becomes smaller than t q if t\ or t 3 is 
examined for smaller time scale than r a . 

In order to obtain the average lifetime of the DH, we 
define the volume of the heterogeneities as 



/>oo />oo 

A hctcl - (k,t 2 ) = dt 3 dhAF^k,^,^,^). 
Jo Jo 



(19) 
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FIG. 9. (a) Waiting time t-z dependence of the integrated 
three-time correlation function Ah e tero(fc, t2)/ Ai leteTO (k, 0) 
with k = 2-7T for various temperatures. The waiting times 
are normalized by r a for each temperature. The solid curve 
is determined by a fitting with the stretched-exponential form 
for each temperature, (b) Average lifetime of DH Thetcro nor- 
malized by the Q-relaxation r a versus temperature T. Inset: 
Relation between two time scales, r^ctoro and r Q . The straight 
line with the slope 1.5 is presented as a viewing guide. 



We examined the t 2 dependence of Ah e tero(fc, h)- Fig- 
ure [Sfa) shows Ahetero(fc)i2)/Ahctoro(fc70) as a function 
of the waiting time t 2 normalized by r Q at each tem- 
perature. From Fig. [9ta), we see that Ahetero rapidly 
decays to zero at higher temperatures and that the time 
scale is comparable to r a . In contrast, at lower temper- 
atures, the relaxation of Ahctcro occurs on a time scale 
larger than r a . Ahctcro (M2)/ Ahctcro (&, 0) can be fitted 
by the stretched-exponential function exp[— (£2 /Thetero) ] > 
where Thetero can be regarded as the average lifetime of 
the DH. We obtain the approximate relation as Thetero — 

/ / Thctcro^l, £3)^1*3/ / / dtidt 3 . We plot Thetero at 

each temperature T in Fig. E[b) . It is found in Fig. HJb) 
that Thetero becomes much larger than t q as the tem- 
perature T decreases. In practice, the lifetime Thetero is 
approximately 6t q ~ 2400 with c ~ 0.5 at the lowest 
temperature of T = 0.289. Furthermore, as seen in the 
inset of Fig. [^b), we observe the strong deviation be- 
tween the two time scales Thetero and r a that follows the 
power low, Thetero ~ Tq, 15 . Similarly, we determined the 
average lifetimes T^° ero and T 1 ™ tero for the mobile and 
immobile parts. We confirm that t 2 dependences of the 
mobile and immobile parts are close to that of the total 
function and that T™° tCTO and T h ™ tero are comparable to 
Thetero for each temperature (data not shown). 



D. 2D spectra of three-time correlation functions 

Through the use of the analogy to the 2D IR spec- 
troscopy, it is of interest to examine the three-time corre- 
lation function AF^k, t\,t 2 , t 3 ) in the frequency domain. 
The Fourier transformed 2D spectrum is obtained by 

AF 4 (k, oji, t 2 ,uJa) 

pOO poo 

= / dt 3 / dt 1 AF 4 (k,t 1 ,t 2 ,t 3 )e lulltl+lu > 3t3 . (20) 
Jo Jo 




10" 5 10 4 10" 3 ^o' 2 10" 1 10" 5 icr 4 10 3 10" 2 10" 1 10" 5 10" 4 10" 3 io 2 10' 1 
a, 




W 5 10~ 4 10~ 3 10~ 2 10"' 10~ 5 10~ 4 10~ 3 10~ 2 10"' 10~ 5 10~ 4 10~ 3 10~ 2 1(T 1 



03; 

FIG. 10. Imaginary parts of 2D spectra of the three-time 
correlation functions; (a) total Q[AF^(k,U)i,t3,u}3)], (b) the 
mobile part S^AF™ ^, cji, £2, W3)], and (c) the immobile part 
^[AFl rn (k,uj 1 ,t 2 ,uj3)] at T = 0.289. Waiting times are varied 
at ti =0, T a , and 3r a from left to right. The wave vector k 
is chosen as k — 2n. The profile is normalized by the peak 
value at ti = for each panel. 



Similarly, the mobile part AF™° and immobile part 
AF™ can be represented in the frequency domain with 
respect to t\ and t 3 . In Fig. [TU1 we present the imagi- 
nary parts of 2D spectra, the total S[AF4(fc, wi, t 2 , uj 3 )] 
(a), the mobile part ^s[AF™°(k,uji,t 2 ,uj 3 )] (b), and the 
immobile part 3[AFj m (fc, wi, t 2 , w 3 )] (c) at T = 0.289 for 
several waiting times. It is seen in Fig. [TUta) that the 
peak of 3[AFi] appears near the detected slower time 

scale, (wi,W 3 ) - (Tr7etero> r hetero)> which is longer-lived 

for larger waiting times. The peak is diagonally elon- 
gated at t 2 = because of the strong correlations be- 
tween two the frequencies uj\ and u 3 . The elongation 
that is directed toward the high frequency side tends 
to be lost because of the loss of the frequency correla- 
tions at larger t 2 . Moreover, the off-diagonal cross peaks 
become pronounced at (wi,w 3 ) ~ (% NG p> T hetero) and 

( T hetero' T NNGp)- Tne time scale t nngp corresponds to 
the peak position of AF^{k, ti, t 2 , t 3 ) at larger t 2 , as seen 
in Fig. EJa). 

Similar behaviors are observed in the 2D spectra of 
the mobile and immobile parts, Q^AF™ ] and ^[AF™]. 
The mobile part 3[AFf 10 ] is rather horizontally elon- 
gated until uj\ ~ T^p because of the coupling between 
^3 — r iietero and ^he higher frequency uii . which is corre- 
lated to the 2D representation in the time domain that 
is seen in Fig. [HJb). Furthermore, the immobile part 
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S[AF| m ] tends to be almost symmetric for the diagonal 
line u>\ = W3, although the 2D profile in time domain is 
largely asymmetric as seen in Fig. |6jc). 



V. CONCLUSIONS AND FINAL REMARKS 

We have investigated the four-point, three-time den- 
sity correlation function to quantitatively characterize 
the temporal structures of the DH. The correlations de- 
tected by the three-time correlation function can be di- 
vided into two parts, mobile and immobile contributions 
determined from the single-particle displacement during 
the first time interval. These 2D representations in both 
the time and frequency domains that are presented over 
a wide range of time scales enable us to explore the cou- 
plings of particle motions. It is shown that the peak po- 
sitions of the mobile and immobile parts are correlated 
to the dominant time scales of the non-Gaussian param- 
eters. These extracted results are not obtainable from a 
one-time correlation function. 

Furthermore, the progressive changes in the waiting 
time allow us to obtain detailed information regarding 
the correlations of motions decay with the time. The 
waiting time dependence of the multi-time correlation 
function shows the existence of the correlations on larger 
time scales between immobile particles. The multi-time 
correlations allow for the quantification of the average 
lifetime of the DH, Thetero, in glass- forming liquids. Our 
analysis can be regarded as an analogue of the multi- 
dimensional nonlinear spectroscopic analysis applied to 
liquids and biological systems to understand ultrafast dy- 
namics, e.g., the transition from inhomogeneous to homo- 
geneous broadening and the couplings between molecular 
motions. 

We have found that the Thetero becomes much slower 
than the a-relaxation time T a when the system is highly 
supercooled. This is due to the long-lived DH at lower 
temperatures. Our findings show that the presence of the 
new time scale Thetero exceeds that of the a-relaxation 
time, T a . These findings are correlated with recent nu- 
merical studies! 19 i 28 ' 56 ~ 59 ' 95 Such strong deviations and 
decouplings between Thotcro and t q when approaching 
the glass transition temperature have been observed in 
some experiments) 96 ' 97 in which the lifetime of a sub- 
ensemble is measured after selective excitation. Recent 
single-molecule experiments that detect the local mobil- 
ity of probe molecules dispersed in glassy materials have 
other relevance to our simulations2£r-i2i. In these ex- 
periments, the lifetime of the dynamical heterogeneity is 
evaluated from the exchange time between mobile and 
immobile regions, which is found to be much slower than 
the structural relaxation time r a near the glass transition 
temperature. 

It is of great importance to examine the relation be- 
tween the length and time scales of the DH in order to 
characterize the relevant spatiotemporal structures. So 
far, various relations such as r ~ f as seen in critical 



phenomenai^i^iS 7 .^ oit~ exp((£/fc B T) c ) based on 
the random first order transition^ have been proposed. 
In these studies, the four-point correlation function is 
used to extract the length scale £ of the DH, where the 
fluctuation in the two-point correlation function 5F{k 1 1) 
is considered as an order parameter as seen in Eq. ([1]). 
On the contrary, the time scale is typically chosen as the 
relaxation time of the two-point density correlation func- 
tion, r a . Here, wc show that the time scale t associated 
with £ should not be r Q but should instead be the average 
lifetime of the order parameter, i.e., Thotcro- This Thetero is 
hidden in the two-point correlation function and unveiled 
when we apply the the four-point, three-time correlation 
function. 

It is worth mentioning that several recent attempts 
have utilized the multi-time correlation function to de- 
tect heterogeneous dynamics. A third-order nonlinear 
susceptibility is theoretically applied to the glassy sys- 
tems Recently this concept has been tested exper- 
imentally to measure the DH.iSi Furthermore, as previ- 
ously mentioned, recent 2D optical techniques^ 9 - - — can in 
principle provide information on heterogeneous dynam- 
ics via 2D representations of multi-time correlation func- 
tions. We hope that our numerical results will be directly 
compared with the experimental ones in the future. 
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